Міністерство освіти і науки України
Національний університет “Львівська Політехніка”
Звіт
до лабораторної роботи №1
з курсу “Чисельні методи математичної фізики”
“ Застосування методу прямих
для розв’язування одновимірних параболічних задач ”
Варіант №19
Львів 2012
Мета роботи: набути практичних навичок під час розв’язування рівнянь у частинних похідних параболічного типу методом прямих, який дозволяє звести такі рівняння до задачі Коші для системи звичайних диференціальних рівнянь з тридіагональною матрицею Якобі.
Постановка задачі
Розв’язати методом прямих задачу:
,
,
,
звівши її до задачі Коші для системи звичайних диференціальних рівнянь з три діагональною матрицею. Для чисельного розв’язування такої системи використати неявні лінійні багатокрокові методи, реалізованими в програмі STIFF.
Введемо рівномірну сітку на інтервалі :
. Замінимо точне значення розв’язку в точці наближеним: , а частинні похідні за змінною різницевими похідними з другим порядком апроксимації, отримаємо систему звичайних диференціальних рівнянь:
.
Крайові умови матимуть вигляд:
.
Використовуючи це значення крайових умов, виключимо їх з першого та останнього рівняння системи, дістанемо задачу Коші для системи звичайних диференціальних рівнянь:
,
Отриману задачу розв’яжемо за допомогою програми STIFF з точністю в режимі MF=25.
Текст головної програми та підпрограм DIFFUN, PEDERV
program test2
implicit real*8(a-h,o-z)
dimension y(50,13),ymax(50),error(50),pw(200),
fsave(100),iwork(50)
common/stcom1/t,h,hmin,hmax,eps,n,mf,kflag,jstart,maxord
common/stcom2/hused,nqused
common/stcom3/ml,mu
common/stcom4/nstep,nfun,njac
common/pcom/coef1,coef2,nm1
nydim=50
eps=1.d-2
kb=3
401 continue
n=50
t=0.0d0
tend=10.d0
h=1.d-3
do 10 i=1,n
10 y(i,1)=0.d0
hmax=tend
hmin=1.d-15
dx=1.d0/dble(n)
coef1=1.d0/dx**2
coef2=1.d0/(2.d0*dx)
nm1=n-1
jstart=0
mf=25
ml=1
mu=1
maxord=5
write(0,20) mf,eps
20 format(//3x,'mf=',i2/,' eps='d11.3)
nstep=0
nfun=0
njac=0
do 30 i=1,n
30 ymax(i)=dmax1(dabs(y(i,1)),1.d0)
40 continue
do 45 i=1,n
45 ymax(i)=dmax1(dabs(y(i,1)),ymax(i))
call stiff(y,ymax,error,pw,fsave,iwork,nydim)
if(kflag.eq.0)go to 60
write(0,50) kflag
50 format(/' kflag=',i2/)
stop
60 continue
if(dabs(tend-t).le.1.d-15) go to 90
if(tend-t-h) 80,40,40
80 e=tend-t
s=e/h
do 85 i=1,n
do 85 j=1,jstart
85 y(i,1)=y(i,1)+y(i,j+1)*s**j
t=t+e
go to 60
90 continue
write(0,556) h,t,(y(i,1),i=1,n)
556 format(1x,5d16.8)
write(0,95) nstep,nfun,njac
95 format(/' nstep=',i4,' nfun= ',i5,' njac=',i4)
kb=kb+1
if(kb.ge.3) go to 402
eps=eps*1.d-2
go to 401
402 continue
stop
end
subroutine diffun (n,t,y,ydot)
implicit real*8 (a-h,o-z)
dimension y(1),ydot(1)
common/pcom/ coef1,coef2,nm1
ydot(1)=coef1*(1.d0-2.d0*y(1)+y(2))-
(coef2*(y(2)-1.d0))**2
ydot(n)=2.d0*coef1*(y(nm1)-y(n))
do 10 i=2,nm1
ydot(i)=coef1*(y(i+1)-2.d0*y(i)+y(i-1))
-(coef2*(y(i+1)-y(i-1)))**2
10 continue
return
end
subroutine pederv(n,t,y,pw,nydim)
implicit real*8 (a-h,o-z)
dimension y(1),pw(1)
return
end Отримані результати
/
Висновок: виконавши цю лабораторну роботу, я навчився розв’язувати рівняння у частинних похідних параболічного типу методом прямих, який дозволяє звести такі рівняння до задачі Коші для системи звичайних диференціальних рівнянь з тридіагональною матрицею Якобі. Для чисельного розв’язування системи звичайних диференціальних рівнянь, я використав лінійні багатокрокові методи, які реалізовані в програмі STIFF в режимі MF=25 (методи Гіра із стрічковою структурою ма...
